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Abstract 


Modifications that allow for more accurate 
modeling of flow fields when strong shocks are 
present have been made to three dimensional 
transonic small disturbance (TSD ) potential 
theory. In addition, Engqulst-Osher type- 
dependent differencing was Incorporated into the 
solution algorithm. The modified theory was 
Implemented and tested In the XTRAN3S computer 
code. Steady flows over a rectangular wing with 
a constant NACA 0012 airfoil section and an 
aspect ratio of 12 were calculated for 

freestream Mach numbers (M«) of 0.82, 0.84, 
and 0.86. Comparisons between results obtained 
using the modified and unmodified TSD theories 
are presented along with results from a three 
dimensional. Euler code. Nonunique solutions In 
three dimensions are shown to appear for the 
rectangular wing as aspect ratio Increases. 
Steady and unsteady results are shown for the 
RAE tall plane model at M« = 0.90. Comparisons 
are made between calculations using unmodified 
theory, modified theory and between experimental 
data. 


Nomenclature 


a* critical speed of sound 

A Mo/ 

AR wing aspect ratio 

B 2M</ 

b semi span 

Cl wing lift coefficient 

C p pressure coefficient 

Cp S entropy generated pressure 

coefficient 


Cp unsteady pressure coefficient 

first harmonic of Cp 
c r reference chord 

c v constant volume specific heat 

E 1 - M. 2 

f(x,t) Instantaneous wing surface 

position 

fj 1*0, 1,2, 3 flux terms In equations 1-5 

F -1/2 (y+ 1)M/ (NASA Ames) 

-1/2 [3 - (2-Y)Mo/jM/ (NLR) 


G 


H 


k 

Mo. 

R 

t 

T 


U. 

u 


a 

Y 

A 

n 


a) 


-1/2 ( y - 3 ) M =. 2 (NASA Ames) 

-1/2 (NLR) 

- (y-1)M. 2 (NASA Ames) 

- M. 2 (NLR) 

reduced frequency (uc r /2U„) 
free stream Mach number 
nondlmensional Ized critical speed 
of sound 

nondlmenslonallzed time in units 
of chords of travel 
ratio of normal velocity down- 
stream of shock to normal 
velocity upstream of shock 
free stream velocity 
streamwlse perturbation velocity 
streanwise, spanwlse, and ver- 
tical distance nondlmenslonallzed 
by reference chord 
angle of attack 
ratio of specific heats 
jump in quantity 
fraction of semispan 
small disturbance velocity 
potential 

angular -frequency, 2wFrequency 


Subscripts 


i 

L 

U 

tip 


free stream conditions 
section lift 
lower surface 
upper surface 
wing tip 


Superscripts 

* sonic value 


Introduction 


When shock waves appear in flow fields, 
aerodynamic loads predicted using potential flow 
theory can be grossly inaccurate or even 
multivalued. The problem of nonuniqueness of 
full potential solutions was first discussed by 
Stelnhoff and Jameson 1 . They found numerical 
evidence that it may be a property of the 
partial differential equation that describes the 
flow field. In two dimensions, potential flow 
theory can produce multiple (nonunique) steady 
state solutions for a fixed angle of attack and 
free stream Mach number. Three steady flows can 
exist, each satisfying the flow equations and 
the associated boundary conditions. This 
phenomenon occurs only when strong embedded 
shocks are present. Further research on the- 
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nonuniqueness of full potential solutions was 
done by Salas and Gumbert^ who demonstrated that 
the problem was not confined to a particular 
airfoil or flow condition. Later, Williams, 
Bland, and Edwards J showed that two dimensional 
transonic small disturbance theory also produces 
multiple solutions. This characteristic Is a 
result of the assumptions of Isentroplc and 
Irrotational flow made when reducing the Euler 
equations to the potential flow equations. 
There have been no reported cases of multiple 
solutions obtained when solving the Euler 
equations. 

Because of the above problems, any 

structural analysis carried out using loads 
predicted using TSD theory cannot be considered 
reliable. The Inaccuracy In the TSD theory is 
due to the way In which the computed shock waves 
are modeled. Classical TSD theory does not 
model the jump In entropy that a fluid particle 
experiences as It passes through a shock wave, 
and calculated shocks can have the wrong 
strength and the wrong location. 

Fuglsang and Williams'* addressed the 

problem of multiple solutions in 2-D TSD theory 
by making changes to the flow equations. They 
did this by modifying the streamwlse flux (to 
Improve the shock jump conditions), modeling the 
effects of the Increased entropy by calculating 
the entropy jump through shocks and convecting 
it downstream, and corrected the pressure 
coefficient to take Into account the effects of 
entropy. The modifications were Implemented and 
tested In the XTRANZL 3 code. The resulting 
method gave results that were In much better 
agreement with Euler solutions than calculations 
made with the unmodified TSD theory. 

The purpose of this paper is to extend 
these nonlsentroplc modifications to three 
dimensional TSD theory. They have been 

Implemented and tested In the XTRAN3S finite 
difference computer code b . The incorporation of 
these modifications Into XTRAN3S, as will be 
shown, permits a more realistic modeling of 
transonic flow fields. 


Governing Equations 

The flow equation currently solved In 
XTRAN3S Is 


3f 0 + ^_ + 3f 2+ j** =0 

at ax ay az 


(l) 


where 


f 0 - -(A* t + B<t> x ) (2) 

^1 = E *x + F *x 2 + ^y 2 (3) 


The constants A, B, and E are defined as 


A 



B = 2Mj 


( 6 ) 


E * 1 - M^ 2 * 

The constants F, G, and H depend on the 
assumptions made In deriving the flow equation. 
In the current effort, two separate sets were 
used. The first set, known as the NASA/Ames 
coefficients 7 . Is 

F « -1/2 (y+DM,, 2 

G - 1/2 (Y-3)M a 2 (7) 

H = - (y-I)M., 2 

These coefficients are obtained when the 
unsteady TSD equation Is derived In a manner 
similar to the steady three dimensional equation 
developed by Lomax et al. 7 . The second set, or 
NLR coefficients/ are 

F « -1/213 - (2-y)M 0 , 2 JH (b 2 
G - - 1/2 M^ 2 (8) 

H = - M ^ 2 


These coefficients were obtained by deriving the 
unsteady TSD equation from the mass conservation 
equation 8 . 


The boundary conditions imposed on the 
field are 

flow 

9=0 

far upstream 

( 9) 

V*t - 0 

far downstream 

(10) 

* z * o 

far above & below wing 

(ID 

* y * 0 

far spanwise & at wing root 

(12) 

a<> 2 = o 

along wake 

(13) 

ACp » 0 

along wake 

(14) 

The condition 
wing surface 

that the flow be tangent to 
is 

the 

*2 - f X + V 

o < X < 1, y < y t1p , z = 0 

(15) 


where + indicates the upper surface, and - 
denotes the lower surface. The Instantaneous 
wing surface position is given by f. 


f 2 = *y + H Vy 


(4) Finite difference solutions of (1) are 

obtained using an alternating direction implicit 

(ADI) method. The ADI method is time accurate 

W and requires one sweep in each of the coordinate 
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directions to advance the solution one time 
increment. In XTRAN3S Version 1 .5* Murman-Cole 3 
type-dependent differencing is used for mixed 
flows. For reasons to be discussed later, 
Engqui st-Osher iU differencing is implemented for 
the present study. 

TSD Modifications 


upstream velocity normal to the shock. 

For lifting cases there will be a jump in 
the upper and lower surface values of C ps . 
Assuming that this jump in C ps , AC ps , is 
convected downstream at free stream speed, the 
following expression is obtained 


This section describes the modifications 
that were made to include noni sentropic effects 
in TSD theory. These changes were first 
implemented for two-dimensional flows*. In the 
present study, they have been extended to three 
dimensions by applying them at each spanwise 
plane as though the flow in those planes was two 
dimensional. This technique still allows for 
the calculation of swept shocks. Details of the 
modifications are given below, although the 
method that is used is the same as that 
presented in Ref. 4. 

Modified Flux 


The original streamwise flux, given by (3), 
has been replaced by 


a ac 
I t 


PS 



( 22 ) 


Equation (22) is integrated to obtain AC ps 
along the wake. The condition of continuous 
pressure across the wake then becomes 


3A$ 3 Aift _ 1 

at ax '2 iL ps 

Without the noni sentropic modifications, 
C ps = 0, and the usual condition 

_ n 

at ax 


(23) 


(24) 


is obtained. 


= (Y + 1 )m£r (V*V - V 2 /2) + G4> y 2 (16) 


where 


R 


V 

V 


r Wl + 


2(1 - Q 

(Y + 1)M 2 


i + U7H + R) 


v(u ) = 


R 2 - 1 
2R 


(17) 


(18) 

(19) 


This new flux term is identical to that used in 
Ref. 4 except for the term which is 

included to model spanwise flow effects. This 
replacement of (3) by (16) requires only a 
change in the streamwise sweep of the ADI 
method. 


Entropy Corrections 

When a fluid particle passes through a 
shock, its entropy increases abruptly and then 
remains constant unless another shock is 
encountered. The Increase In entropy results in 
a pressure decrease proportional to the strength 
of the shock. This nondimensional pressure 
loss, C ps , is given by 

- 

C = — 

PS Yfr-DM. 2 (20) 


where the entropy jump is found using the 
Rankine-Hugoniot relationship. 



where T is the ratio of the downstream to 


Finally, the entropy generated pressure is 
added to the linear isentropic term 


C = C + C 

p ’Linear ps 


(25) 


where 


C 


p Linear 


-2(* x + V 


(26) 


Engquist Osher Switching 


While testing the noni sentropic 
modifications it was found that solutions would 
sometimes converge with velocity overshoots 
immediately upstream of shock waves. Because 
the resulting spike in the pressure distribution 
occured only at the shock, the differencing 
method was changed. Using Engqui st-Osher 
differencing eliminated the velocity overshoots 
upstream of shocks. Although not demonstrated 
in three dimensions, Engqui st-Osher differencing 
has been shown to dissipate expansion shocks 
and to improve code stability^. 


Results and Discussions 


Flows past an infinite aspect ratio wing 
with a NACA 0012 section were calculated to 
determine if the 3-D method yields the 2-0 

solution for that case. For all such cases 

involving the NACA 0012 section the NLR 
coefficients were used. The effect of aspect 

ratio on the solution for a rectangular wing 
using a NACA 0012 airfoil section is presented. 
A wing of aspect ratio 12 was choosen for a 
parametric study - varying Mach number and angle 
of attack - to determine the effects of the 
noni sentropic modifications. Comparisons with 
solutions of the Euler equations are shown. 

Comparisons of steady and unsteady results for 
the RAE tailplane model with experimental data 
are also shown. 
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NACA 0012 


.5 


2-0 Results 

The two dimensional XTRAN3S results were 
obtained by calculating flow past an infinite 
aspect ratio rectangular wing. This study was 
performed In order to verify that using XTRAN3S 
In this manner produces the same multiple 
solutions as XTRAN2L.* 4 

Figure 1 shows the pressure distribution on 
a NACA 0012 airfoil at M« 3 0.84 and angle of 
attack a = o°. This Mach number and angle of , 
attack were chosen because these conditions are 
in the middle of the 2-0 nonuniqueness region 
for this airfoil. The three dimensional code 
is seen to produce the 2-D nonunique solution. 

The solution has converged to steady state after 
being perturbed by a Gaussian pulse in angle of 
attack. The differences In the upper and lower 
shock locations cause the large lift associated 
with this case. 

The region of Mach numbers in which 
multiple solutions have been found to occur for 
the NACA 0012 airfoil is shown in Fig. 2. Two 
sets of results are shown, those from XTRAN2L 
and the 2-D results of XTRAN3S. Comparisons 
between these two curves show that for the NACA 
0012 airfoil the nonuniqueness regions predicted 
by the two codes are very similar, Differences 
in the two results can be attributed to 
differences in computational grids, differenc- 
ing methods, and far field boundary conditions. 

Figure 3 shows lift-coefficient time 
histories that result from a Gaussian pulse in a 
of one-fourth degree at M* = 0.84. The angle 
of attack returns to zero by 150 chords of 
travel. The curve calculated using the 
unmodified code reaches a steady state value of 
0.45. Replacing the original streamwise flux, 

Eq. (3), with the modified flux of Eq. (16) 



x/c 


Fig. 1 Steady state pressure distribution for 

NACA 0012, M - .84, a * 0° using XTRAN3S 
simulating 2-D flow. 



Fig. 2 Steady state 2-D lift values for the 
NACA 0012 at a = 0°. 



Fig. 3 Lift time histories for the NACA 0012 
using a pitch pulse, amplitude = 1/4°, 

M = .84 with XTRAN3S simulating 2-D flow. 


causes the lift to converge to a value of 0.3. 
When the entropy corrections, Eqs. (22) - (25), 
are included with the modified flux, the lift 
returns to zero after the pulse in angle of 
attack. These results are similar to those 

obtained by Fuglsang and Williams. 4 

3-D Results 

To examine the effects of finite aspect 
ratio on the flow solution, calculations were 
made for flow conditions at which multiple 
solutions exist in two dimensions. For this 
study a rectangular, AR = 12 wing with a NACA 
0012 airfoil section was used. Throughout this 
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discussion the above wing will be referred to as 
"Wing A". 

Figure 4 shows lift histories for various 
aspect ratios using the same Gaussian pulse 
described earlier at = 0,84. These results 
are from the unmodified XTRAN3S code. Increas- 
ing the aspect ratio leads to nonzero lift- 
coefficients similar to those found in two 
dimensions, though at lower aspect ratios 
the lift returns to zero. Thus, nonunique 
solutions of the TSD potential equation can 
occur in three dimensional flow. 

Figure 5 shows the steady state lift 
coefficient versus the Inverse of the aspect 
ratio. The results for the rectangular wing 
with a NACA 0012 section at * 0.84 and 
a * 0° show that multiple solutions exist for 
aspect rati os greater than about 20 for the 
unmodified code. Though wings of such* high 
aspect ratios are rare, these results 
demonstrate that 3-D potential flows have the 
same problems with multiple solutions as 2-D 
flows If the wing aspect ratio is large enough. 



Fig. 4 Lift time histories for the NACA 0012 
using a pitch pulse on various aspect 
ratio wings. 


The problem concerning the velocity 
overshoots when using Murman-Cole differencing 
is shown in Fig. 6. The cause of the overshoot 
is unknown although it was not encountered until 
the flux modifications were added to the XTRAN3S 
code. The overshoot occurs at various locations 
along the span. The use of Engquist-Osher 
differencing eliminated the problem as shown in 
Fig. 7. 

Figures 8-14 illustrate the effect of the 
present modifications for an aspect ratio 12 



Fig. 5 Steady state lift coefficient versus 

aspect ratio for a rectangular wing with 
a NACA 0012 airfoil section at M * .84, 

a = 0°. 


1.25 r 



-.75 h 


-1.25 1 I I I l l 

0 .2 .4 .6 .8 1.0 

x/c 

Fig. 6 Steady state pressure distribution using 
Murman-Cole switching. 


wing with a NACA 0012 section at M w = 0.82, 
0.84, and 0,86. These Mach numbers were chosen 
because of the 2-0 results shown in Fig. 2. 
Figures 8-10 are for M - 0.82, which is outside 
of the 2-0 nonuniqueness region shown in Fig. 
2. Figure 8 shows a comparison of steady Ci 
vs. a calculated using (1) unmodified XTRAN3S, 
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Fig. 7 Steady state pressure distribution using 
Engqulst-Osher switching. 



Fig. 8 Total lift coefficient versus angle of 
attack for wing A at M = .82. 


(2) modified XTRAN3S, and (3) the FL057MG Euler 

code. 2 * * * * * * * * 11 In this case, the modified XTRAN3S 

results agree well with the Euler calculations 
at a = 1°. The lift curve slope calculated with 
the modified XTRAN3S code had a value 
approximately half that calculated using the 
unmodified XTRAN3S. Figure 9 shows a comparison 
of spanwise lift distribution at 1° angle of 
attack. The modified TSO theory shows excellent 
agreement with the Euler results except at the 
tip where the Euler code seems to have some 
difficulty, believed to be due to gridding near 
the tip. The corresponding comparison of the 
pressure distributions is shown in Fig. 10 for 

1° angle of attack. Only two lower surface 

pressure distributions are shown since the 

effects of the strong shock modifications 

decrease toward the tip in the same manner as on 

the upper surface. The excellent agreement in 

spanwise load distribution is a result of the 

more accurate shock location. 


Figures 11-13 show correspondi ng results 
for^ M«, = 0.84 which is in the middle of the 
nonuniqueness region shown in Fig. 2. At this 

higher Mach number, as shown in Fig. n, the 
lift curve calculated with the original code has 
a very steep slope initially and then levels off 
(due to the shock reaching the trailing edge). 
The nonisentropic modifications produce a nearly 
linear lift curve, as does the Euler code - 
although the agreement with the Euler code is 
not as good as at M. = 0.82. 

Figures 12 and 13 show section lift 

distributions at a = 1° and 2° for M„ = 0.84. 
The present modifications greatly improve the 
agreement between TSD results and Euler results 
in a way similar to that shown In Fig. 9. 
Figure 13 also shows how the section lift from 
the unmodified code begins to diverge from the 
elliptical behavior exhibited by the Euler 
results. The divergence in total lift 

coefficient for the nonisentropic and the Euler 
results as the angle of attack increases is 
believed to be due to assumptions of TSO theory 
(thin wing, small angle of attack) which are 
less valid as the angle of attack Increases. 

When the Mach number is increased to 
0.86 unmodified TSD theory predicts the upper 
surface shock location to be at the trailing 
edge for any significant angle of attack. 
Figure 14 shows that the total lift curve using 
the modified code is much more linear and closer 
to the Euler results although the 1'lft values 
are high. 

RAE 

The nonisentropic theory was tested on the 
aspect ratio (AR = 2.41), tapered swept wing of 
the RAE tailplane model. Theairfoil used is 
approximately the NACA 64A010.2. In order to 
compare with results from Ref. 13, the NASA Ames 
coefficients were used along with the grid 
transformation used by Bennett 1J . 



Fig. 9 Section lift coefficients versus 

semispan for wing A at H = .82, a = 1°. 
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Fig. 13 Section lift coefficients versus 

semi span for wing A at M « .84, a =■ 2°. 



Fig. 14 Total lift coefficient versus angle of 
attack for wing A at M « .86. 


Steady results for M - 0.90, a = -.3° with 
and without the nonlsentroplc modifications are 
shown in Fig. 15. As the tip of the wing is 
approached, the shock Increases in strength and 
the nonlsentroplc modifications have their 
largest effect. Because the angle of attack Is 
small for this case, the entropy convection has 
little effect whereas the new flux causes the 
largest change in the results of Fig. 15. 

Unsteady results for the above case are 
shown in Fig. 16. For this case, the wing was 
oscillated at 70 Hz about a = -.3° with an 
amplitude of .57°. The reduced frequency used 
was k * .44. The center of rotation was 68 

percent chord. As these results show the 
noni sentroplc modifications result in a slight 
improvement in agreement with experimental 


data. The effects are greatest near the tip. 
The effects of entropy convection is greater 
here because the amplitude of oscillations is 
greater and creates a difference in shock 
strengths on the upper and lower surfaces. 


Conclusion 

Nonlsentroplc modifications to three 
dimensional transonic small disturbance (TSD) 
theory have been formulated. The modifications 
were incorporated and tested in the 3-D TSD 
code, XTRAN3S, for flows past a rectangular NACA 
(1012 wing of aspect ratio 12. A comparison of 
Euler solutions and the results obtained using 
the modified theory shows that a significant 
improvement was made in modeling the Euler 
solutions. One of the most noticeable 
improvements Is the more accurate lift curve 
slope. Also discovered is the property of the 
unmodified code to produce multiple solutions 
for wings of sufficiently large aspect ratios. 

The nohlsentropic modifications were also 
tested on RAE tall plane model for steady and 
unsteady flows. Comparisons of modified, 
unmodified, and experimental results showed that 
the modified results Improved the agreement with 
experiment for both steady and unsteady results 
particularly near the tip where the shock is the 
strongest. 
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